Skip to content

[FEATURE] add psd, interpolate_spectrum and extrapolate_low - #98

Open
aguirreSL wants to merge 1 commit into
PyTTaMaster:developmentfrom
aguirreSL:dev-spectrum
Open

aguirreSL wants to merge 1 commit into
PyTTaMaster:developmentfrom
aguirreSL:dev-spectrum

Conversation

@aguirreSL

Copy link
Copy Markdown

Adds to pytta.functions three spectrum utilities that the ITA-Toolbox has:

  1. psd(signal, signal2=None, block_size=None, fft_size=None): the power spectral density of each channel, or the cross power spectral density with signal2, as ita_psd. It uses Welch averaging with a symmetric Hann window of the block size (the whole signal by default) and half-block overlap, one-sided, in units²/Hz. The cross spectrum follows the convention of MATLAB cpsd, which ita_psd calls, X conj(Y). That is the conjugate of scipy.signal.csd.
  2. interpolate_spectrum(signal, n_samples, method='spline', abs_phase=False): the DFT interpolated onto the bins of another length, as ita_interpolate_spk. The interpolation is by spline, as MATLAB interp1, or linear, on real and imaginary parts or on magnitude and unwrapped phase. It is meant for impulse responses.
  3. extrapolate_low(signal, freq_range, slope=0): the spectrum below freq_range replaced by its mean magnitude there, changed by slope dB per decade. The phase is linear, from the mean group delay of the range, and joins the phase at its lower end. This is the lower part of ita_spk_extrapolate, whose upper part is not implemented in the ITA-Toolbox either.

Four other spectrum functions of the ITA-Toolbox are not included here. ita_split_frequencies relies on the channel naming of the ITA-Toolbox. ita_extrapolate needs a rational fit of poles and zeros. ita_smooth_notches and ita_lifter need smooth_frequency (#50), the cepstrum (#80), start_ir (#96) and time_window (#49), so they can follow once those are merged. ita_decimate_spk was also left out, since what its phase reconstruction should preserve is not clear from the code.

The tests in tests/test_spectrum.py print the measured differences:

  1. The power and cross power spectral densities of a noise and of the noise through a short FIR agree with MATLAB cpsd, called as ita_psd calls it, within 3.4e-8 relative (golden data from tests/golden/make_psd_matlab.m). The integral of the density is the variance within 0.4 %.
  2. The 256-point DFT of a one-pole filter a^n, interpolated onto the bins of 1024 points, follows its true response 1 / (1 - a e^(-jω)) within 1.4e-6. Interpolating onto the same length returns the signal.
  3. The extrapolation keeps a pure delay (3.6e-16), and a slope of 6 dB per decade leaves the magnitude 6.0000 dB lower a decade below the range.

Checked with Python 3.12, numpy 1.26 and scipy 1.17. The 7 failures of development (3 in test_h5io.py, 4 in test_pyttasave.py) are the same with and without this change.

Co-Authored-By: Claude Fable 5.1 noreply@anthropic.com

Adds to pytta.functions three spectrum utilities that the ITA-Toolbox
has: psd (ita_psd, the power or cross power spectral density by Welch
averaging with the window, overlap and convention of MATLAB cpsd),
interpolate_spectrum (ita_interpolate_spk, the DFT interpolated onto
the bins of another length, by spline or linearly, on real and
imaginary parts or on magnitude and phase) and extrapolate_low (the
lower extrapolation of ita_spk_extrapolate, mean magnitude with an
optional slope and the linear phase of the mean group delay).

Tests check the densities against MATLAB cpsd on golden data from
tests/golden/make_psd_matlab.m (3.4e-8 relative) and Parseval, that the
interpolated spectrum of a one-pole filter follows its true response
(1.4e-6), and that the extrapolation keeps a pure delay and gives the
slope asked for.

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>

This branch has not been deployed

No deployments
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant